function Collect_Fox_discrete_FEonly_grID_wtp_tomlab_varagin_j_bs(subsample,inc,set_seed,matlab_seed,get_nb_draws)

%matlab_seed=1353;
rng(matlab_seed)

	tmp_file1=strcat('/Users/shoude/Dropbox/eegap/EEgap_data_code_heter_SM/Data/Matlab_estimation/choiceset_identifier_trimester_week_zipcode_2008_2012_v11022017_struct_',num2str(subsample),'_',num2str(inc),'_seed_',num2str(set_seed),'.csv');
    id=load(tmp_file1); 
	Nmax=size(id,1);
	if subsample == 11000
		bs_max=10*ceil(Nmax/10000);
    elseif subsample == 44000
		bs_max=10*ceil(Nmax/20000); 
    end

%Temporary code
bs_max=16;

c=1
for j=1:bs_max	

	try	
		results_file=strcat('/Users/shoude/Dropbox/eegap/EEgap_data_code_heter_SM/Results_Estimation/demand_est_grID/','Fox_discrete_mixedDemoshrt_FEonly_WTP_v',num2str(subsample),'_inc',num2str(inc),'_sd',num2str(set_seed),'_cty_tau1_grid_v3_','bs_',num2str(j),'_rep.mat');
	    fox=load(results_file);
	    j 
	    eta=fox.param_v(:,1);
	    theta=fox.param_v(:,2);
	    tau(:,c)=fox.param_v(:,3);
	    weight(:,c)=fox.beta;
	    c=c+1;
	catch
		c=c;	
	end
    
end

if subsample==11000
	tau_v={'tau1','tau2','tau3','tau4','tau5','tau6','tau7','tau8','tau9','tau10','tau11'};
	weight_v={'weight1','weight2','weight3','weight4','weight5','weight6','weight7','weight8','weight9','weight10','weight11'},;
elseif subsample==44000
	tau_v={'tau1','tau2','tau3','tau4','tau5','tau6','tau7','tau8','tau9','tau10','tau11','tau12','tau13','tau14','tau15','tau16'};
	weight_v={'weight1','weight2','weight3','weight4','weight5','weight6','weight7','weight8','weight9','weight10','weight11','weight12','weight13','weight14','weight15','weight16'},;
end

if subsample==11000
	colNames = {'eta','theta','tau1','tau2','tau3','tau4','tau5','tau6','tau7','tau8','tau9','tau10','weight1','weight2','weight3','weight4','weight5','weight6','weight7','weight8','weight9','weight10'};
elseif subsample==44000
	colNames = {'eta','theta','tau1','tau2','tau3','tau4','tau5','tau6','tau7','tau8','tau9','tau10','tau11','tau12','tau13','tau14','tau15','tau16','weight1','weight2','weight3','weight4','weight5','weight6','weight7','weight8','weight9','weight10','weight11','weight12','weight13','weight14','weight15','weight16'};
    %colNames = {'eta','theta','tau1','tau2','tau3','tau4','tau5','tau6','tau7','tau8','tau9','tau10','tau11','tau12','weight1','weight2','weight3','weight4','weight5','weight6','weight7','weight8','weight9','weight10','weight11','weight12'};
end

%colNames = {'eta','theta','tau1','tau2','tau3','tau4','tau5','tau6','tau7','weight1','weight2','weight3','weight4','weight5','weight6','weight7'};
Table_beta = [eta,theta,tau,weight];
Table_beta_n = array2table(Table_beta,'VariableNames',colNames);
			
results_all=strcat('/Users/shoude/Dropbox/eegap/EEgap_data_code_heter_SM/Results_Estimation/demand_est_grID/','Fox_discrete_mixedDemoshrt_FEonly_WTP_v',num2str(subsample),'_inc',num2str(inc),'_sd',num2str(set_seed),'_cty_tau1_grid_v3_all_bs','.csv');
    
writetable(Table_beta_n,results_all);

%Create an average across bs

dim_weight = size(weight);
pos_weight = zeros(dim_weight);
weight_id  = find(weight>1e-3);
pos_weight(weight_id)=1;

avg_weight = sum(weight,2)./sum(pos_weight,2);
inf_id = isinf(avg_weight);
avg_weight(inf_id) = 0;
Table_avg_weight = [eta,theta,avg_weight];
%Table_avg_weight = array2table(avg_weight); 
results_all=strcat('/Users/shoude/Dropbox/eegap/EEgap_data_code_heter_SM/Results_Estimation/demand_est_grID/','Fox_discrete_mixedDemoshrt_FEonly_WTP_v',num2str(subsample),'_inc',num2str(inc),'_sd',num2str(set_seed),'_cty_tau1_grid_v3_avg_weight','.csv');
   
dlmwrite(results_all,Table_avg_weight,'precision','%10.6f');

				
end